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This paper presents a comparison between two imple- 
mentations of the Ffowcs Williams and Hawkings equation 
for airframe noise applications. Airframe systems are gen- 
erally moving at constant speed and not rotating, so these 
conditions are used in the current investigation. Efficient 
and easily implemented forms of the equations applica- 
ble to subsonic, rectilinear motion of all acoustic sources 
are used. The assumptions allow the derivation of a sim- 
ple form of the equations in the frequency-domain, and the 
time-domain method uses the restrictions on the motion to 
reduce the work required to find the emission time. The 
comparison between the frequency domain method and the 
retarded time formulation reveals some of the advantages 
of the different approaches. Both methods are still capable 
of predicting the far-field noise from nonlinear near-field 
flow quantities. Because of the large input data sets and 
potentially large numbers of observer positions of inter- 
est in three-dimensional problems, both codes utilize the 
message passing interface to divide the problem among 
different processors. Example problems are used to demon- 
strate the usefulness and efficiency of the two schemes. 

Nomenclature 

c 0 ambient speed of sound 

/ integration surface defined by / = 0 

F-i dipole source terms 

H Heaviside function 

i \f^l 

Mi local source Mach number vector, V{/c 
M Mach number, \Mj\ 

hi Outward directed unit normal vector 
p pressure 

Q monopole source term 

Qi components of vector in Eq.(6) 

Qn Qi'ft'i 

ri radiation vector, Xj — & 

r magnitude of radiation vector, |r. t | 

t time 

Ui Cartesian fluid velocity components 
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Vi Cartesian surface velocity components 
x , y, z Cartesian observer coordinates 

Greek: 

, a2 mapped surface coordinates 
p x/1 - M 2 

5 (f) Dirac delta function 

5 jj Kronecker delta 

p fluid density 

r retarded or emission time, t — r/c Q 

£, 77, C, Source coordinates 

Superscript: 

perturbation quantity (e.g. p' = p - p a ) 
unit vector 
time derivative 

Subscript: 

e emission or retarded quantity 

o freestream quantity 

r inner product with f; 

r et quantity evaluated at retarded time r 

„ inner product with hi 

Introduction 

Despite recent advances in computational aeroacous- 
tics, numerical simulations that resolve wave propagation 
from near-field sources to far-field observers are still pro- 
hibitively expensive and often infeasible. Integral tech- 
niques that can predict the far-field signal based solely on 
near-field input are a means to overcome this difficulty. 

The Ffowcs Williams and Hawkings 1 (FW-H) equation 
is a rearrangement of the exact continuity and Navier- 
Stokes equations. The time histories of all the flow vari- 
ables are needed, but no spatial derivatives are explicitly 
required. The solution to the FW-H equation requires a 
surface and a volume integral, but the solution is often well 
approximated by the surface integral alone. Singer et al? 
and Brentner and Farassat 3 have shown that when the sur- 
face is in the near field of a solid body, the FW-H approach 
correctly filters out the part of the solution that does not 
radiate as sound, whereas the Kirchhoff method produces 
erroneous results. Many applications of the FW-H and 
Kirchhoff methods can be found in the area of rotorcraft 
acoustics. 4-7 The FW-H method has typically been applied 
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by having the integration surface coincide with solid bod- 
ies, but the method is still applicable when the surface is 
off the body and permeable. The codes developed in this 
work are valid for both cases. 

For three-dimensional flows, the time-domain FW-H 
formulations developed by Farassat 8 are efficient and 
amenable to numerical computations. Some simplifica- 
tions are applied here based on restrictions of the surface 
motion, but the development follows that of Farassat very 
closely. The frequency-domain version of the FW-H uses 
the form of the equations developed by Lockard 9 for two- 
dimensional problems. However, the derivation was done 
using Cartesian tensor notation, and is equally valid in 
three-dimensions as long as the appropriate Green function 
is employed. 

The remainder of the paper describes the time- and 
frequency-domain formulations in enough detail to show 
the similarities and differences between the methods. A 
discussion of various parallelization strategies follows. Fi- 
nally, two example problems are used to demonstrate the 
utility and efficiency of the schemes. The last example 
involves a noise calculation based on a large-scale CFD 
computation of a landing gear. 


Governing Equations 

The FW-H equation can be written in differential form 10 


as 


9^_ 2 9 2 

9t 2 C ° dxidxi 

9 2 


H(f)ft') = 


dxidxj 

9 


TijH(f) 


- o^[mf)) + ^(QS(f) 


9 


9t 


( 1 ) 


where 


F; 


Pij (MiU.j “f Pij C o p $ij , 


df 


9xj ’ 


Q = ( PoVi + p(Ui - Vi) j . 


( 2 ) 


Pij + pui(uj - Vj) ) — , and (3) 


(4) 


The contribution of the Lighthill stress tensor, Ty, to the 
right-hand side is known as the quadrupole term. The 
dipole term involves an unsteady force, and Q gives rise 
to a monopole-type contribution that can be thought of as 
an unsteady mass addition. The function / = 0 defines 
the surface outside of which the solution is desired. The 
normalization |V/| = 1 is used for /. The total density 
and pressure are given by p and p, respectively. The fluid 
velocities are u,, while the v t represent the velocities of 
the surface /. The Kronecker delta, <Y,j, is unity for i = j 
and zero otherwise. The ambient speed of sound is denoted 


by c 0 . A prime is used to denote a perturbation quantity 
relative to the free-stream conditions denoted by the sub- 
script o. The Cartesian coordinates and time are x.j and 
t, respectively. The usual convention, which is followed 
here, involves a quiescent ambient state with / prescribed 
as a function of time so that it always surrounds a moving 
source region of interest. H(f) is the Heaviside function 
which is unity for / > 0 and zero for / < 0. The deriva- 
tive of the Heaviside function H'(f ') = 6(f) is the Dirac 
delta function, which is zero for / ^ 0, but yields a finite 
value when integrated over a region including / = 0. The 
inviscid part, P, j = p5/j, of the compressive stress tensor 
Pij is used in this work. 

Equation 1 is typically solved using a Green function 
technique. The temporal and spatial convolution of the 
free-space Green function with the source terms yields the 
solution for p’ . In the farfield, p' = <? 0 p' . The three- 
dimensional Green function 10 is 


G(Xi,t-,£i,T ) 


5{t - t+r/c 0 ) 
4nr 


( 5 ) 


where = Xi — & is the radiation vector, and r = r. ( |. 
Because of the delta function, the solution to equation 1 
can be manipulated into various forms, some of which are 
amenable to numerical solution such as formulation 1A of 
Farassat. 8 An alternative approach is to solve the FW-H 
equation in the frequency domain. In this work we fol- 
low the frequency-domain approach of Lockard 9 that is 
restricted to uniform, rectilinear motion, but other methods 
can be used that still allow general motion. In the following 
sections, the current implementations of these methods are 
discussed. 


Time-Domain Method 

A concise, time-domain solution to equation 1 that is ap- 
plicable to nondeforming, porous surfaces can be obtained 
from the derivation of Farassat 11 using the variables 


Qi = ( po ~ p)vi + pui, and 

Li = phi + pui(uj — Vj)hj (6) 

as proposed by Di Francescantonio. 6 The integral solution 
is given by 


/ 

/= o 

+ 


4irp'(xi,t) = 

Qjhj Qjhj(rM r + c a (M r — M 2 )) 


r( 1 - M r ) 2 

/b 

/= o 


r 2 (l - M r 


( 7 ) 

ds 


Ljf-j 


+ 


Lj v j Li A Li 


r(l — M r ) 2 r 2 (l — M r ) 2 


ds 


+ 


/ 

/= o 


LjTj(rMjfj + c 0 (MjTj — M 2 )) 

c 0 r 2 (l - M r ) 3 


ds + p' Q 
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where j/q is the quadrupole term, which is neglected in this 
work. Mi = Vi/c a is the local Mach number vector, and the 
superscripted caret ( A ) denotes a unit vector. The outward 
directed unit normal vector is given by hi. The subscript 
ret denotes evaluation at the retarded or emission time r. 
For uniform, rectilinear motion, the normals are not func- 
tions of time, so finding Q n = Qihi and Qihi is sufficient 
because d/dt(Q n ) = hid/dt(Qi) + Qid/dt(hi) = Qihi. 
The dot above a variable indicates a time derivative. Hence, 
only four variables and their time derivatives are required 
to compute the kernel functions in equation 7. The time 
derivatives are computed using fourth-order, central differ- 
ences. In the current application, the coordinates of the 
surface and the flow field data are known at the grid points 
of meshes on a series of surface patches that comprise the 
complete, closed surface surrounding all sources. The nor- 
mals are determined using the grid metrics on the patches, 
which are computed using fourth-order, central differences. 
Because each surface patch has two independent indices, 
only two variables (ai,a 2 ) are required to describe the 
variation on the mapped surface. The normal is computed 
from the outer product 

n, = {dx/dcti,dy/dai,dz/dai} x 

{dx/da 2 ,dy/da 2 ,dz/da 2 }. (8) 

Whether the normal points inward or outward is dependent 
on the ordering of the data. The current code reads a file 
which either specifies a source point within each surface 
or the correct sign of the normal on each surface. When 
the surfaces come from CFD data, it is usually possible to 
deduce which direction the normal will point. However, a 
visual inspection of the normals is always a good practice 
because incorrect signs are common and can cause errors 
that are difficult to recognize. 

One of the more computationally intensive parts of a 
time-domain calculation is the determination of the re- 
tarded time t. For general motion, a root-finding technique 
must be employed. However, for uniform, rectilinear, sub- 
sonic motion, the Garrick triangle 12 uniquely determines 
the emission time as 

r e = 1 _ r M2 \ M cos (0) + \J 1 - M 2 sin 2 (0)^ 

t = t- r e /c 0 ( 9 ) 

where r e is the distance between the source and the ob- 
server at the time of emission t. The angle between the 
surface velocity vector v r and the radiation vector r* is 6. 
The cos {6) can be determined from the inner product of the 
vectors. The variables Q„ , and their time derivatives 
must then be interpolated to the retarded time. This is an 
intensive part of the calculation because the interpolation is 
performed for every grid point, retarded time, and observer 
position. Changing the indexing of data to increase cache 


reuse, loop unrolling, and inlining of the interpolation rou- 
tine decreased the total run time by half. 

Once the kernel functions in equation 7 have been deter- 
mined, Gaussian quadrature is used to perform the spatial 
integration over the surface. Because this integration is 
performed repeatedly for each observer and retarded time, 
parts of the integral that are independent of time and ob- 
server position are precalculated and stored. The integra- 
tion is performed by linearly interpolating the data from 
the four corners of each cell using three-dimensional shape 
functions commonly employed in finite elements. 13 Defin- 
ing the mapped coordinates as (oq, a 2 ), the elemental area, 
dS, is the magnitude of the normal as specified in equa- 
tion 8. The value of the elemental area at each Gaussian 
quadrature point is then stored for reuse. The values of the 
shape functions at each of the quadrature points are also 
precalculated and stored. Typically, only four quadrature 
points are required over each cell. Four points will integrate 
quadratic functions exactly. Although the function being 
integrated may be of higher order, a quadratic approxima- 
tion is usually sufficient provided the surface shape is well 
behaved. Additional quadrature points and higher-order re- 
constructions within each cell have not been found to be 
very beneficial. Part of the reason may be the oscillatory 
nature of the kernel function. The integration of a linear ap- 
proximation of a sinusoid is often not much different from 
the quadrature of a higher-order fit. Grid refinement is typi- 
cally much more influential on the solution quality than the 
order of the integration. 

Because of the need to sample the source data at the re- 
tarded time, one cannot arbitrarily choose observer times. 
A separate code was written to determine the valid range of 
observer times based on a given range of input data. The 
Garrick triangle 12 can be used to uniquely determine the 
minimum and maximum reception times using equations 9 
to solve for t given r. 

Impenetrable Surfaces 

In many practical applications, the data for the FW- 
H solver is obtained on solid surfaces, and the equations 
can be greatly simplified thereby decreasing the computa- 
tional and memory requirements. The simplifications are 
especially advantageous for the uniform, rectilinear motion 
case. For impenetrable surface data, Uj = Vj and equations 
6 simplify to 

Qi = PoVi and Li = phi . (10) 

Note that Q is independent of time. Furthermore, the nor- 
mals hj are only a function of space. Hence, only the time 
history of the pressure and its time derivative are needed 
when the viscous terms are neglected in Pij. All of the 
source terms are linear, yet the method has been shown to 
give correct results when the surface data is in the nonlin- 
ear near field. 2 The failings of the Kirchhoff method for 
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problems with solid surfaces is clearly not due to its source 
terms being linear, but rather the assumptions used to de- 
rive the equation. 

Frequency-Domain Method 

Several frequency- domain formulations of the FW-H 
equation have been reported in the literature. 14 Here, the 
method proposed by Lockard 9 for two-dimensional flows 
will be extended to three-dimensions. Although this de- 
velopment is restricted to surfaces in rectilinear motion at 
constant speed, it is still useful for airframe noise, where 
these assumptions are usually valid. The motion of the sur- 
face is assumed to be governed by / = /(x + Ut) where 
the components of U are constant velocities describing the 
motion of the surface. An application of the Galilean trans- 
formation from (x, t) to (y, t), 

yi = Xi + Uit, t = t, 

d _ d d _ d d 

dxi dyi dt dt + % ’ 

to equation 1 leads to the convected wave equation 


+ 


d 2 

dy,dy i 


d 2 

r=s + u i u i 


a 2 


dt 




( 12 ) 


2Uj 


a 2 


a 2 

dyidt C ° dyi dy { 


H(f)p ') = 


a 


dyi 


TijH(f) l Fi5 (f) +-[Q n S(f) 


a 


dt 


equation 12 becomes 


^7 ^F i (y,uj)5(f) S j - iujQ n (y,uj)5(f) 

(15) 


The wavenumber is defined by k = w/c 0 , the Mach num- 
ber M = U/c a , and complex number i = \f^ T. Note that 
the transform has been applied to the groupings Ty , F , t , and 
Q because the equation is linear in these terms. However, 
the desirable properties of the FW-H are maintained be- 
cause all of the nonlinear products are included before the 
transformation is applied. In a numerical implementation, 
the products are formed first, and then a fast Fourier trans- 
form (FFT) is applied. As a caution, the FFT must use the 
sign convention of equations 14, or the derivation must be 
modified appropriately. The Green function for equation 
15 when M < 1 can be obtained from a Prandtl-Glauert 
transformation. Denoting the three-dimensional source co- 
ordinates as (£,V,0, and the observer position as ( x,y,z ) the 
Green function 15 is 


G(x,y,z-t,r 1 ,Q=— d 


exp (-ik(d-Mx)/f) 2 ) 


where, after the transformation, the F, and Q n become 

Fi — Fi ) ( Uj + Uj) + PoUiU^Jnj, 

Qn + Ui) ( 13 ) 

The Lighthill stress tensor Ty is unchanged, and / = /( y) 
is now only a function of the spatial coordinates. The sur- 
face velocities u* have been replaced by —Uj, which can 
be inferred from inspection of /(x + Ut) = 0. Note 
that this implies that the mean flow is in the positive direc- 
tion (or equivalently that the surface moves in the negative 
direction) when Ui > 0. The Vj used throughout the time- 
domain formulation are of opposite sign and should not be 
confused with the Uj. Equation 12 is now in a convenient 
form to perform the Fourier analysis. With application of 
the Fourier transform pair 

oo 

F{q(t)} = q{uS) = J q(t ) exp(—iojt)dt and 

— OO 

OO 

•7 7-1 {?(“>)} = q{t) = j" q(u)exp(iwt)du, (14) 


where 

x = (x — £) cos a cos + (y — r]) sin a 
+ (z — C) cos a sin cf>, 
y = —(x — tl)sinacos0 + (y — ri)cosa 

+ (z — () sin a sin <j>, 
z = — (x — £) sin cf> + (z — 0 cos <fr, 

and 

d = \Jx 2 + /3 2 (y 2 + z 2 ) (16) 

The angles are defined such that taxup = W/U, sin a = 
V/M, and M = yj U 2 + F 2 + W>/c 0 . The Prandtl- 
Glauert factor is /? = y/1 — M 2 . The solution to equation 
15 for M < 1 can now be written as 

H(f)c 2 0 p'( y,w) = - J Fi(£,uj) dG ^’® ds 
f = o 

- J *wQ„(£,w)G(y;£) ds 

~ / (17) 

f>0 
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As has already been stated, the quadrupole contribution, 
represented by the last term in equation 17, is typically ne- 
glected because its influence is often small. Furthermore, 
the calculation is somewhat involved and expensive. Cer- 
tain flows exist where the quadrupole cannot be ignored, 
such as those containing significant refraction of waves by 
shear layers and wakes. As long as the integration sur- 
face is placed outside of all regions where T tJ is large, 
the quadrupole contribution is substantially included by the 
surface sources even though the quadrupole integration is 
not performed. This is also true for the time-domain for- 
mulation. 

The frequency-domain solution process involves calcu- 
lating the surface normals and forming the products in F 
and Q for all time at each point on the surface just as in 
the time-domain version. However, the F> and Q n func- 
tions are Fourier transformed, and the surface integrations 
are performed for each frequency of interest instead of for 
each observer time. An inverse FFT can be used to recover 
the acoustic signal in the time domain. For truly periodic 
problems one merely uses a single period of the flow data as 
input to the FW-H code. However, for more complicated, 
aperiodic flows, windowing the data is required. The win- 
dowing should be applied to F and Q n after their mean 
values are subtracted. The subtraction has no effect on the 
calculated noise because the derivatives of G all contain 
uj, and equation 17 shows that there is no contribution to 
the noise at u = 0 when the quadrupole tenn is neglected. 
The minimal amount of time data typically available from a 
computational aeroacoustics calculation may lead to some 
inaccuracies in the windowed FFT, but short time records 
are often just as much of an impediment for time-domain 
formulations because information about the frequency con- 
tent is usually desired. 

When the input to FW-H code is from a harmonic, lin- 
earized Euler solver, F and Q n should also be linearized 
because the amplitudes from the linearized code may not 
be physical. If they are too large, the nonlinearities in 
the source terms can produce erroneous results. One must 
be careful when performing the linearization because the 
perturbation velocities u,- are not necessarily small. For in- 
stance, on a solid surface u; = — (/, . Only a minor change 
in the code allows one to have a single code that is use- 
ful for input from linear and nonlinear flow solvers. The 
frequency domain approach is particularly efficient for har- 
monic data because only a single frequency needs to be 
calculated, and the FFT’s do not need to be performed. 

A disadvantage of this particular frequency-domain for- 
mulation is that the source and observer are always a fixed 
distance apart, and all Doppler effects are lost. In a time- 
domain calculation, the distance between the observer and 
the source can be changed for each time step to simulate 
a flyover condition. Most CFD computations and exper- 
iments are carried out in a laboratory frame with the ob- 


server distances fixed, so this is not a major issue. How- 
ever, comparisons with actual flight data should include the 
Doppler effects. 

Impenetrable Surfaces 

As was shown for the time-domain formulation, the 
frequency-domain version can also be significantly simpli- 
fied when the input data is obtained on solid surfaces. For 
impenetrable surface data, u j = — Uj and equations 13 
simplify to 

Q = -poUihi and F, = phi. (18) 

Note that Q is steady in time and has no impact on the 
frequency domain solution. Hence, only the time history 
of the pressure is needed. One only needs to determine 
the Fourier transform of the pressure and scale that result 
by the appropriate normal to obtain the Ft terms. The 
savings in memory and computational requirements are so 
great that the solid surface formulation should be employed 
whenever possible. In the current implementation, differ- 
ent versions of two subroutines are called depending on 
the case being run. Although some additional coding is re- 
quired, run times can be reduced by 60% and the memory 
load by over 70%. 

Parallel Implementation 

Although one normally thinks of acoustic analogy com- 
putations as being extremely efficient, calculations involv- 
ing long time records and many observers can quickly be- 
come expensive. The cost of computing the time history 
at a single point using the FW-H may actually exceed that 
of a standard CFD method. However, the FW-H approach 
allows the observer location to be anywhere outside of the 
source region, whereas the CFD method must have grid 
points from the source to the observer. The cost involved 
with all those intermediate points and the errors incurred 
in the long range propagation make standard CFD an in- 
appropriate choice for most far-field noise computations. 
Still, when mapping out the directivity in three-dimensions, 
hundreds of observer locations may be required. In the re- 
alistic problem of a landing gear given in the examples sec- 
tion, the FW-H computation of a single observer using the 
porous surface formulation requires seven minutes on an 
SGI 250MHz R10K processor. The input record contains 
4096 time steps at 82,219 grid points, which consumes 1.4 
GB of disk space. Mapping out a directivity is a time con- 
suming process when performed serially. This motivated 
the modification of the code to use the Message Passing 
Interface (MPI) to perform parallel computations. Other 
researchers 16 have used MPI in conjunction with acoustic 
analogy methods. 

FW-H solvers are ideal candidates for parallelization be- 
cause the calculation of the signal at each observer is in- 
dependent, and the contributions from each portion of the 
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data surface combine linearly. Even the computations at 
each time step or frequency are independent. Hence, one 
has many choices of how to split the problem. Some initial 
testing was done with different processors dealing with dif- 
ferent observers. However, this requires that each processor 
have access to the entire data record for all the surface 
patches. In the landing gear case, 181 patches comprise 
the total surface. Either all of the nodes have to read all 
the data, or one node must read it and broadcast it to all 
the others. The total record for the landing gear is expected 
to be over 10 GB on a medium mesh, and over 60 GB for 
a fine mesh. Either reading or passing that much data is 
not reasonable. This same problem occurs when different 
processors handle different frequencies or time steps. The 
other option is to divide the problem by surface patches, 
and sum all of the contributions at the end. Each proces- 
sor only needs to read the data for the particular patches 
assigned to it. Hence, the data is only read once and only 
passed if the data is not directly accessible by all the nodes. 

Two paradigms were used to investigate the parallel im- 
plementation. First, a standard load-balancing approach 
was used. The size of each of the patches was read, and 
the largest patch assigned to the processor with the least 
points until all the patches were assigned. This strategy is 
commonly employed in multiblock CFD codes. This ap- 
proach worked well when all of the processors were identi- 
cal and dedicated. However, because the code only passes 
a very minimal amount of data, it is typically run over a 
standard network on the second processor of SGI Octane 
workstations in the lab. These machines vary in speed, 
and sometimes both processors are in use when the job 
starts. Occasionally, one of the processors would take much 
longer than all of the others to complete. To circumvent 
this problem, the master-slave paradigm employed by Long 
and Brentner 17 was used. Here, one node does no work 
other than to tell requesting nodes what patch they should 
process. Again, the largest patches are assigned first so 
that small patches are being used when the job nears com- 
pletion and the variability in processor speed is important. 
This paradigm worked extremely well and has resulted in a 
nearly linear speed up on an SGI Origin. Although the mas- 
ter node does not do any useful work, it needs to respond 
quickly to the requests from worker nodes to keep them 
from becoming idle. The master’s job can be assigned to 
a slow node or to a dual processor machine that also has a 
slave process. 

In a heterogeneous environment where one is trying to 
use idle machines, the master-slave paradigm is preferred. 
Furthermore, the load-balanced approach and master-slave 
approach only affect a single subroutine in the code, so it is 
very easy to switch between the two depending on the local 
operating environment. 


Test Examples 

Monopole in Flow 

As a first demonstration of the current implementations 
of the three-dimensional FW-H equation, the field from 
a monopole source is computed in the far field using the 
present technique. The source moves in the —x direction at 
Mach 0.5. An equivalent flow involves a fixed source at the 
origin in a uniform flow in the +x direction. The complex 
potential for the monopole is given by Dowling and Ffowcs 
Williams 15 as 

4>(x, y , z, t) = A-^~- exp >(^-k(d-Mx)/0 2 ) (19) 

The variables needed in the FW-H equation are obtained 
from the real parts of p' = —p 0 (d(j>/dt + U 0 d(j>/dx), «• = 
d<j>/dxi, and p' = p'/ c o- Equation 19 is written in a 
laboratory frame where the flow is moving over a station- 
ary source. The source terms in the FW-H equation are 
calculated from the flow variables evaluated over two pe- 
riods on the surface. For this case, M = U 0 /c 0 = 0.5, 
ujI/co = 47t/46, A/(lc 0 ) = 0.01 and the integration sur- 
face is a cube that extends from —51 to 51 in all three coor- 
dinate directions. The reference length is l. Fifty uniformly 
spaced points are used on each side of the box. Figure 1(a) 
compares the directivity from the calculations to the ana- 
lytic solution in the 2 = 0 plane. Figure 1(b) makes a 
similar comparison for the time history at (50/, 0, 0). The 
agreement is excellent, demonstrating that the formulations 
are valid for problems with a uniform mean flow. Simi- 
lar agreement was found when each of the six sides of the 
cube comprising the data surface were deformed to look 
like a Gaussian bell provided enough points were provided 
to adequately resolve the variation. The calculations were 
performed in single precision, which appears sufficient as 
long as the acoustic signal on the surface is not less than 
five orders of magnitude smaller than the mean. A lack 
of smoothness in a time-domain formulation solution is an 
indication of a precision problem. 

Landing Gear 

This example involves the calculation of the noise gener- 
ated by the unsteady flowfield surrounding a landing gear. 
An aerodynamic and acoustic analysis of a similar land- 
ing gear was performed by Souliez et a/. 18 The near field 
in this problem is highly nonlinear, and different parts 
shed vorticity at different frequencies. Figure 2(a) shows 
an instantaneous snapshot of the perturbation pressure on 
all of the solid surfaces. The freestream Mach number 
is 0.2, and the gear is mounted on a flat plate. The ref- 
erence length is the diameter of the wheels which is 3.7 
in (0.09398 m). The input data for the acoustic calcula- 
tion is obtained from a three-dimensional, time-dependent 
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(a) Directivity at r = 50 1 



(b) Time history at (50 1, 0, 0) 

Fig. 1 Solution comparisons for a monopole in a M = 0.5 
flow. The pressure is nondimensionaiized by p 0 Co 

CFD solution using the code CFL3D. 19 - 20 CFL3D was de- 
veloped at NASA Langley Research Center to solve the 
three-dimensional, time-dependent, thin-layer Reynolds- 
averaged Navier-Stokes (RANS) equations using a finite- 
volume formulation. The CFD data used in this work is 
discussed in more detail in the paper by Li et al. 2X 

The noise calculation involves 181 total patches com- 
prising the data surface. All of the patches are on solid 
surfaces, so only the pressure histories are needed. 147 
patches are on the gear itself and 31 are on the plate above 
the gear. So far, over 24,000 nondimensional time samples 
of l/c 0 At. = 0.02 have been collected, but only half of the 
data is sufficiently free of transients to be useful for acous- 
tic calculations. The transient problem is exacerbated by 


high-frequency sources that are growing through most of 
the time record. Some of the waves from these sources can 
be seen on the door in Figure 2(a). Figure 2(b) shows the 
time history of the pressure on the oleo in the contraction 
just below the door as indicated by the arrow in Figure 2(a). 
The time history shows the complex, intermittent charac- 
ter of the signal. The large amplitude oscillations occur 
around one kHz (model scale) at varying intervals and ir- 
regular durations. Riding on top of the signal are high 
frequency oscillations generated by resonances in small, 
triangular shaped spaces between the yokes and the door. 
These spaces are found on the upper and lower yokes in 
both the upstream and downstream junctures with the door. 
Unfortunately, the signals took a long time to saturate, so 
the calculation had to be run much longer than anticipated 
to eliminate the transients. The upstream and downstream 
cavities are slightly different in size, so the resonances oc- 
cur at 20.5 kHz in the upstream cavity and at 25.4 kHz in 
the downstream one. 

Beyond the intermittency in the signal, another compli- 
cation for obtaining spectral information about the solution 
is the predominantly low frequency content of the gear. 
Very long sampling times are needed to resolve these fre- 
quencies. Although the time-domain code can exactly cal- 
culate the signal at the observer, the lack of long sampling 
times is a problem when the spectral content is needed. 
Nonetheless, the comparison of the spectra at an observer 
directly beneath the gear from the two FW-H codes is en- 
couraging. The observer is 121 away from the gear. Figure 
3(a) shows that the spectral content is nearly identical. In 
both cases, a hamming window was applied to minimize 
the effects of the aperiodic signal on the FFT’s. The time 
histories in Figure 3(b) are also similar. The curves are 
offset because the frequency-domain calculation does not 
include the effect of the zero or steady mode. All of the 
primary features of the signal can be observed in both re- 
sults. Some discrepancy should be expected because the 
effect of the window is already included in the frequency 
domain results because the window is applied to the input 
data. However, the window is applied to the time-domain 
results after the calculation, so its effect is not seen in Fig- 
ure 3(b). 

Three different data records were used to investigate the 
influence of the duration and the particular time interval 
used. Because of the complexity of the time histories on 
the gear, some variation should be expected, but drastic 
changes would indicate an improper sample length or that 
the flow has yet to eliminate transients. These sorts of vari- 
ations were observed in many of the initial calculations. 
However, the comparison in Figure 4 shows only minor dif- 
ferences between the three calculations. A time sample of 
8192 was used, which represents the time required for the 
flow to pass by a wheel 32 times. The full record was used, 
then it was subdivided into two records of 4096 samples. 
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(a) Perturbation pressure on solid surfaces 


(a) Spectra 
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(b) Time history on the oleo 

Fig. 2 CFD results for a landing gear in a M = 0.2 flow. 

The comparison between the results for the full record and 
the latter record is very good. Somewhat more discrepancy 
is observed with the results for the first 4096 samples, but 
it is still within the variation that should be expected for 
such a complex flow. Even though there are not as many 
samples available from the computation as one can obtain 
from an experiment, the power spectral density should be 
calculated by averaging over the different samples to obtain 
a single answer, which should be stationary. 

One of the concerns with the CFD calculation was the 
accuracy of the solution on the wall above the landing gear. 
Although the plate doesn’t produce any significant fluid 
mechanic fluctuations that would act as sound sources, it 
does reflect acoustic signals, and some of the vorticity shed 
by the side bars interacts with the wall. The grid on the 


(b) Time history 

Fig. 3 FW-H results for an observer directly below a landing 
gear in a M = 0.2 flow. 

plate in the CFD calculation coarsens very quickly away 
from the gear, so the accuracy of high frequency signals is 
likely to be very poor over much of the plate. To investigate 
the impact of the accuracy on the plate, the FW-H solvers 
were run with the wall ignored and the gear surfaces mir- 
rored about the plate. Figure 5 shows the comparison with 
results obtained by including and excluding the wall. At 
lower frequencies, the case including the wall has much 
higher levels than either the no wall or mirrored cases in- 
dicating that the signal on the wall is an important source 
of noise. The results in the figure are from the frequency- 
domain code, hut very similar results were obtained from 
the time-domain code. Figure 5(b) zooms in on the high 
frequency tones. For the tones, the wall and no wall re- 
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Fig. 4 FW-H results for an observer directly below a landing 
gear in a M = 0.2 flow. Calculations from frequency-domain 
code comparing effect of sample length. 

suits are nearly identical which supports the idea that the 
high-frequency signals on the wall are artificially damp- 
ened rapidly so that they have little impact on the acoustic 
calculation. The mirrored result is lower at 20.5 kHz and 
significantly higher at 25.4 kHz. The variability in the 
effect of mirroring with frequency is likely caused by in- 
terference effects which would be minimized for observers 
farther from the gear. 

Because of the complexity of the landing-gear geometry 
and the various locations of the sources, one would expect 
the observer location to be important. Figure 6 compares 
the spectra for observers at three locations. All of the ob- 
servers are 12/ from the gear. The first location is directly 
beneath the gear. The other two observers are at the same 
vertical and streamwise locations as the wheels, but on op- 
posite sides. The observer on the oleo side doesn’t have 
the door in the way to obscure some of the sources on the 
oleo and yokes. This is most evident for the 25.4 kHz tone 
in Figure 6(b). Both side observers see much stronger tone 
signals indicating that the yoke resonances primarily radi- 
ate horizontally. Part of the reason that the strength of the 
tones observed underneath the gear is much less may be 
because of blockage by the truck and wheels. The 20.5 
kHz tone has the same amplitude for both side observers. 
The reason for the lack of influence of the door on this tone 
has not been investigated. At lower frequencies, the signal 
beneath the gear has a larger amplitude indicating that the 
wheels and truck are best observed from below. The por- 
tion of the signal around 1 .4 kHz is completely absent for 
the side observers. As expected, there is a significant direc- 
tivity associated with the gear, but the complex inteiplay 
of source location, geometry, and observer position needs 
further investigation. 



(a) 



(b) 

Fig. 5 FW-H results for an observer directly below a landing 
gear in a M =0.2 flow. Results from frequency-domain code 
comparing the effect of including the wall and mirroring the 
data. 

Performance 

In the landing gear calculations, all of the frequencies 
were calculated in the frequency domain code, and an 
equivalent number of time steps in the time domain code. 
Table 1 shows timing results for the serial runs with 4096 
samples and all 181 surfaces. These calculations involved 
solid surface data, so the impenetrable surface simplifica- 
tions were employed. Run times are more than doubled if 
the surface is assumed to be penetrable. All of the frequen- 
cies were calculated by the frequency-domain code, and an 
equivalent number of observer times in the time-domain 
code. When all of the data is computed, the frequency- 
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(a) 



(b) 

Fig. 6 FW-H results for an observer directly below a landing 
gear in a M = 0.2 flow. Results from frequency-domain code 
comparing the effect of observer location. All observers are 
12 / from the gear. The side observers have the same vertical 
location as the wheels. 

domain code is considerably faster. However, the ordering 
of the input data was kept the same for both codes so they 
could share input files. The ordering used is more natural 
for the frequency-domain code which must perform FFT’s 
on the time records. The time-domain code must interpo- 
late the source terms to the retarded time, and this is done 
much more efficiently when the data is ordered such that 
the variable index varies fastest. The current code rear- 
ranges the data internally to obtain the preferred ordering. 
Even accounting for this extra work, the frequency-domain 
code would still be more efficient. Furthermore, the gear 


spectra is primarily concentrated at low frequency and at 
a few select tones. The number of frequencies computed 
could be significantly reduced in the frequency domain 
code. However, fewer temporal samples could be calcu- 
lated in the time-domain code as long as the tones could be 
neglected. The time-domain code also has the advantage 
that one only needs to calculate the new portions of the ob- 
server signal as new input data is made available. Because 
the FFT’s are performed on the input for the frequency- 
domain code, there is a minimum sample length necessary 
to produce reasonable results. 


Code 

CPU Time (s) 

Time Domain 

1089 

Frequency Domain 

425 


Table 1 CPU time comparison of serial computations of a 
landing gear signal for one observer from 4096 time samples 
on 181 patches with 82,219 grid points. Calculations per- 
formed on a 250 MHz, RlOk SGI Octane. 


Code 

Observers 

CPU Time (s) 

Time Domain 

1 

88.0 

3 

128.1 

Frequency Domain 

1 

56.5 

3 

85.0 


Table 2 CPU time comparison of parallel computations of 
a landing gear signal from 4096 time samples on 181 patches 
with 82,219 grid points. Sixteen SGI Octane worker nodes of 
nonuniform speed were used in the master-slave calculations. 

For problems the size of the landing-gear calculation, 
the computations are normally performed in parallel using 
MPI. Because of the high utilization of machines dedi- 
cated for parallel calculations, the master-slave paradigm 
was used to take advantage of the second processor on six- 
teen SGI Octanes in the lab. The processor speeds varied 
from 195 to 250 MHz. The master node was attached to 
a RAID system that stored the input data. Table 2 shows 
the timing results for calculations involving one and three 
observers. The speedup from the serial case is clearly not 
linear, but the run times have been significantly reduced. 
Furthermore, much of the time is spent doing file input 
and output. The master node is not capable of reading the 
data fast enough to satisfy all of the worker nodes. In this 
problem, a total of 1 .4 GB of data must be read and trans- 
ferred across the network. When the code first starts, all 
of the worker nodes simultaneously try to access their data 
from the master. Distributing the data before the calcula- 
tion would be more efficient, but having the data in one 
place is much more convenient. Furthermore, as the num- 
ber of observers increases, the work to disk access ratio 
improves resulting in better efficiency as can be seen from 
the three observer cases in table 2. The point of this exam- 
ple is not to show that perfect efficiency has been obtained, 
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but rather to demonstrate that the MPI implementation has 
enabled directivity mappings of thousands of observer lo- 
cations to be performed relatively quickly on machines that 
are typically idle. 

Conclusions 

Two different formulations of the three-dimensional 
Ffowcs Williams and Hawkings equation for sources in 
uniform, subsonic, rectilinear motion have been presented. 
They are efficient enough to be used to perform far-field 
predictions from large data sets provided by either compu- 
tation or experiment. Comparisons of the solutions for two 
example problems showed excellent agreement between 
the frequency- and time-domain formulations. Although 
the frequency- domain version of the code is somewhat 
faster, both algorithms are efficient enough to be used for 
computations involving large data sets when MPI is used 
to distribute the problem to many computers. Each formu- 
lation has different advantages that make it more attractive 
for certain classes of problems, but our tests have shown 
that they produce basically equivalent results for the prob- 
lems investigated. 
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